
airport_fn <- file.path(wkdir,'local','gis_data','018_transportation','Airports','airports2008globalGIS.csv')
airports_df = read.csv(airport_fn)
roads_afr_fn <- file.path(wkdir,'local','gis_data','018_transportation','roads','rds_remi_afr.shp')

require(Rcpp)
output_version <- '2020-09-21' #Sys.Date()
options(rasterTmpDir='C:/temp/')

##
## Euclidean distance to airport in any country in Africa
## Euclidean distance to airport in same country
## Travel time to airport in same country
## Travel time to airport in any country in Africa
##

# adm2 has admin codes ie GID_0
ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_adm2_x_d01') %>%
  st_set_crs(.,4326)

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,st_drop_geometry(ntl.zones))

# airports 2008
air <- st_as_sf(airports_df,coords=c("longitude","latitude"),crs=4326)
air_join <- st_join(air,subset(ntl.zones,select=c(OBJECTID,GID_0)))
air_join4 <- subset(air_join,GID_0 %in% c("CMR","NER","NGA","TCD"))

ntl.zones.pts["dist2airport_afr"] <- unlist(nngeo::st_nn(ntl.zones.pts, air, k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts["dist2airall"] <- unlist(nngeo::st_nn(ntl.zones.pts, air_join4, k = 1, returnDist = T)$dist) / 1000

## within country
ntl.zones.pts["dist2airport_cmr"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(air_join,GID_0=="CMR"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts["dist2airport_ner"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(air_join,GID_0=="NER"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts["dist2airport_nga"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(air_join,GID_0=="NGA"), k = 1, returnDist = T)$dist) / 1000
ntl.zones.pts["dist2airport_tcd"] <- unlist(nngeo::st_nn(ntl.zones.pts, subset(air_join,GID_0=="TCD"), k = 1, returnDist = T)$dist) / 1000

## add identifiers to grid ##
dist.df <- subset(as.data.frame(ntl.zones.pts),select=c(OBJECTID,GID_0,dist2airport_afr,dist2airall,dist2airport_cmr,dist2airport_ner,dist2airport_nga,dist2airport_tcd))
save.dta13(dist.df,paste0(wkdir,'/proc_data/AIRPORTDIST',output_version,'.dta'))

##
## travel time
##

#Tracks, improved, paved, highways and speed (r) but do Euclidean distance only (d) => two measures will be different if two measures are not strongly correlated. 
#Highway = 80. paved = 60. improved = 40. dirt = 12. no road = 6
#Assumption = free crossing (open border) (o) (later on, we will have p = porous, p1 p2.)
#Sigma = 1, 2, 3, 4 (1 2 3 4)
#For example, baseline = lcmacro4 for each cell-year


## speeds in kph
roadcat2speed <- function(road_cat) {
  
  if (road_cat==8) {
    return(80)
  }
  if (road_cat==1) {
    return(60)
  }
  if (road_cat==3) {
    return(40)
  }
  if (road_cat==5) {
    return(12)
  } 
  if (road_cat==0) {
    return(12)
  } else {
    return(6)
  }
}

fishnet_tif <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_x_d01.tif')
fishnet_grid <- raster(fishnet_tif)

##
## roads
##
road_lk_chad_shp = paste0(wkdir,'/local/gis_data/018_transportation/roads_panel')

if (!file.exists(paste0(road_lk_chad_shp,'.shp'))){
  roadsafr_sf <- st_read(dsn=dirname(roads_afr_fn),layer=basename(roads_afr_fn)) %>% 
    sf::st_make_valid()
  
  reg_sf <- st_read(dirname(adm0_lake_chad_shp), basename(adm0_lake_chad_shp))
  GID_0 <- reg_sf$GID_0
  
  roadsafr_prj_sf <- st_transform(roadsafr_sf,crs(reg_sf))
  roads <- subset(roadsafr_prj_sf, COUNTRY %in% countries)
  
  write_sf(roads,dirname(road_lk_chad_shp),basename(road_lk_chad_shp),driver='ESRI Shapefile')
} else {
  roads <- read_sf(dirname(road_lk_chad_shp),layer=basename(road_lk_chad_shp))
}   

# slope in percent
Slope_globe <- raster(file.path(wkdir,"local/gis_data/006_elevation/meanslope/meanslope"))
crs(Slope_globe)<- CRS("+init=epsg:4326") 
Slope_crop <- crop(Slope_globe, alignExtent(ntl.zones,Slope_globe,snap='out'))
# convert 100
Slope_radians = Slope_crop / 100
Slope <- Slope_radians*pi/180

##
## return ((6 * exp(-3.5 * abs(x + 0.05))) * 0.6)
## waldo in km/h

offpath <- function(x) {
  return ((6 * exp(-3.5 * abs(x*pi/180 + 0.05))) * 0.6)
}

## radians to degrees
onpath <- function(x) {
  return (6 * exp(-3.5 * abs(tan(x) + 0.05) ) )
}

grid_res_km = 11.1

# input: slope in degrees 
# output: km / h 
slope_kph <- raster::calc(Slope, onpath)

# cost in hours per cell
# time = distance / speed
# raster::area no weights
slope_hr.prj <- projectRaster(slope_kph,fishnet_grid, res = 0.1, method='bilinear')
slope_hr.prj <- mask(slope_hr.prj,fishnet_grid)
slope_hr.cell = sqrt(raster::area(slope_hr.prj)) / slope_hr.prj

# by period
roads.t = subset(roads,select=c(COUNTRY,R2008))
# transform road classes into travel cost to cross a pixel
# grid resolution in km
# max speed for each cell 
roads.t$speed = unlist(lapply(roads.t$R2008,FUN=roadcat2speed))

if(!file.exists(roads.grd.fn)){
  roads.grd.t = rasterize(roads.t,s,field = 'speed', fun='max', background=NA)
  writeRaster(roads.grd.t, roads.grd.fn, driver= "GTiff",overwrite=TRUE)
} else {
  roads.grd.t <- raster(roads.grd.fn)
}

# cost in hours per cell
# time = distance / speed
# raster::area no weights
roads.hr.cell = sqrt(raster::area(roads.grd.t)) / roads.grd.t

costgrid.stack <- stack(slope_hr.cell,roads.hr.cell)
# least cost of road path and natural path
cost.grd.t <- stackApply(costgrid.stack, c(1,1), fun='min',na.rm=TRUE) 

## update cost-distance grid with background levels from waldo tobler function
## the values in the input raster are costs. 
## take the arithmetic mean of the cost first and then take the reciprocal of that to get the conductance. 
## ie. The traveller experiences half of the cost of the origin cell and half of the cost of the destination cell, (cost1 + cost2)/2.
cost.t = gdistance::transition(cost.grd.t, transitionFunction=function(x) 1/mean(x), directions=8)

##
## END COST
##


##
year <- 2008
#for (country in countries){
#  print(paste0("country : ",country))

calcPoints <- function(acled_loc){
  # add small value when no ntl #
  acled_loc[acled_loc==0] = 1
  rast.df <- as.data.frame(rasterToPoints(acled_loc))  
}

air_join4 <- subset(air_join,GID_0 %in% c("CMR","NER","NGA","TCD"))
objectdid_airports <- unique(st_drop_geometry(air_join4["OBJECTID"]))
ntl.zones.airports <- subset(ntl.zones, OBJECTID %in% unlist(t(objectdid_airports)))

## ALL 4 countries
tt_cost2airport.hour <- gdistance::accCost(cost.t, as.matrix(st_coordinates(air_join4)))
print(paste0("Travel time (open border): ",names(tt_cost2airport.hour)))
writeRaster(tt_cost2airport.hour,paste0(wkdir,'/proc_data/tt2airport.tif'),driver='GTiff',overwrite=TRUE)
ntl.zones$tto2air <- exact_extract(tt_cost2airport.hour,ntl.zones,fun='mean')

## CMR
air_cmr = subset(air_join,GID_0=="CMR")
tt_cost2aircmr.hour <- gdistance::accCost(cost.t, as.matrix(st_coordinates(air_cmr)))
print(paste0("Travel time (open border): ",names(tt_cost2aircmr.hour)))
writeRaster(tt_cost2aircmr.hour,paste0(wkdir,'/proc_data/tt2aircmr.tif'),driver='GTiff',overwrite=TRUE)
ntl.zones$tto2aircmr <- exact_extract(tt_cost2aircmr.hour,ntl.zones,fun='mean')

## NER
air_ner = subset(air_join,GID_0=="NER")
tt_cost2airner.hour <- gdistance::accCost(cost.t, as.matrix(st_coordinates(air_ner)))
print(paste0("Travel time (open border): ",names(tt_cost2airner.hour)))
writeRaster(tt_cost2airner.hour,paste0(wkdir,'/proc_data/tt2airner.tif'),driver='GTiff',overwrite=TRUE)
ntl.zones$tto2airner <- exact_extract(tt_cost2airner.hour,ntl.zones,fun='mean')

## NGA
air_nga = subset(air_join,GID_0=="NGA")
tt_cost2airnga.hour <- gdistance::accCost(cost.t, as.matrix(st_coordinates(air_nga)))
print(paste0("Travel time (open border): ",names(tt_cost2airnga.hour)))
writeRaster(tt_cost2airnga.hour,paste0(wkdir,'/proc_data/tt2airnga.tif'),driver='GTiff',overwrite=TRUE)
ntl.zones$tto2airnga <- exact_extract(tt_cost2airnga.hour,ntl.zones,fun='mean')

## TCD
air_tcd = subset(air_join,GID_0=="TCD")
tt_cost2airtcd.hour <- gdistance::accCost(cost.t, as.matrix(st_coordinates(air_tcd)))
print(paste0("Travel time (open border): ",names(tt_cost2airtcd.hour)))
tt_cost2airtcd.hour <- mask(tt_cost2airtcd.hour,ntl.zones)
writeRaster(tt_cost2airtcd.hour,paste0(wkdir,'/proc_data/tt2airtcd.tif'),driver='GTiff',overwrite=TRUE)
ntl.zones$tto2airtcd <- exact_extract(tt_cost2airtcd.hour,ntl.zones,fun='mean')


## add identifiers to grid ##
tt.df <- subset(as.data.frame(ntl.zones),select=c(OBJECTID,GID_0,tto2air,tto2aircmr,tto2airner,tto2airnga,tto2airtcd))
save.dta13(tt.df,paste0(wkdir,'/proc_data/AIRPORT',output_version,'.dta'))